function maps = calcPixelMaps(data, varargin)
% Calculate all kinds of orientation and direction pixel maps
%
% inputs
% data: registered data stack
% 'stim_on': indices of stimulated frames either
%            [0 0 0 1 1 1 0 0 0 3 3 3 0 0 0 2 2 2 ..] where 0 indicate the
%            non-stimlated frames and increasing integer numbers the stim
%            conditions or
%            [-1 -1 -1 90 90 90 -1 -1 -1 0 0 0 ..] where -1 indicate the
%            non-stimlated frames and stim conditions are set by value of
%            angle
% 'stim_off': logical indices of baseline frames (optional)
% 'tune_max'
% max values for scaling
% 'ColorShift'
% 'OutputPath'
% 'fileID' : experiment ID

maps = struct;
colorShift = 60/360;
maps.tune_max = 0.4;
outputpath = [];
stim_on = [];
stim_off = [];
stim_off_set = 0;
lowpass = 2; % (um)
usesmoothing = 1;

if ~isempty(varargin)
    numIndex = find(cellfun('isclass', varargin(1:end-1), 'char'));
    for ind = 1:length(numIndex)
        switch lower(varargin{numIndex(ind)})
            case 'stim_on'
                stim_on = varargin{numIndex(ind) + 1};
                if min(stim_on) == -1
                    gap = unique(diff(unique(stim_on(stim_on > -1))));
                    stim_on = fix(stim_on / gap + 1);
                end
                if ~stim_off_set
                    stim_off = (stim_on == 0);
                end
            case 'stim_off'
                stim_off = varargin{numIndex(ind) + 1};
                stim_off_set = 1;
            case 'colorshift'
                colorShift = varargin{numIndex(ind) + 1};
            case 'tune_max'
                maps.tune_max = varargin{numIndex(ind) + 1};
            case 'outputpath'
                outputpath = varargin{numIndex(ind) + 1};
            case 'lowpass'
                lowpass = varargin{numIndex(ind) + 1};
                if lowpass == 0
                    usesmoothing = 0;
                    lowpass = 2;
                end
        end
    end
end

if isempty(stim_on)
    error('Stimulus indices not specified. At least ''stim_on'' has to be set.')
end

if ~exist(outputpath, 'dir') && ~isempty(outputpath)
    mkdir(outputpath);
    mkdir([outputpath filesep 'dir']);
    mkdir([outputpath filesep 'ori']);
end

pixel_in_x = size(data,2);
fov_in_um = 200;
lowpass_pix = lowpass * pixel_in_x / fov_in_um;
kernel_size = ceil(lowpass_pix * 5);
sp_filter = fspecial('gaussian', kernel_size, lowpass_pix);    % spatial lowpass filter in pixel unit (after binning)

[maps.base, maps.base_sm] = CalcBaseImage(data, stim_off, sp_filter);
maps.dir_F = CalcFImages(data, stim_on);
[maps.dir_dF, maps.dir_dF_sm, maps.ori_dF, maps.ori_dF_sm] = CalcdFImages(maps.dir_F, maps.base, sp_filter);

if usesmoothing
    [maps.dir_ratio, maps.ori_ratio] = CalcdF_FImages(maps.dir_dF_sm, maps.ori_dF_sm, maps.base_sm);
    [maps.max_dir_dF, maps.max_dir_ratio_change, maps.max_ori_dF, maps.max_ori_ratio_change] = ...
        write_dF_Images (maps.dir_dF_sm, maps.ori_dF_sm, maps.dir_ratio, maps.ori_ratio, 'outputpath', outputpath);
    maps.dir_dF_params = CalcMapParams (maps.dir_dF_sm);
    maps.dir_ratio_params = CalcMapParams (maps.dir_ratio);
    
    maps.ori_dF_params = CalcMapParams (maps.ori_dF_sm);
    maps.ori_ratio_params = CalcMapParams (maps.ori_ratio);
else
    [maps.dir_ratio, maps.ori_ratio] = CalcdF_FImages(maps.dir_dF, maps.ori_dF, maps.base);
    [maps.max_dir_dF, maps.max_dir_ratio_change, maps.max_ori_dF, maps.max_ori_ratio_change] = ...
        write_dF_Images (maps.dir_dF, maps.ori_dF, maps.dir_ratio, maps.ori_ratio, 'outputpath', outputpath);
    maps.dir_dF_params = CalcMapParams (maps.dir_dF);
    maps.dir_ratio_params = CalcMapParams (maps.dir_ratio);
    
    maps.ori_dF_params = CalcMapParams (maps.ori_dF);
    maps.ori_ratio_params = CalcMapParams (maps.ori_ratio);
end


if ~isempty(outputpath)
    maps.dir_angle = GenerateAngleMap (maps.dir_dF_params.th, 'colorshift', colorShift, 'outputpath', [outputpath '\dir']);
    maps.ori_angle = GenerateAngleMap (maps.ori_dF_params.th, 'colorshift', colorShift, 'outputpath', [outputpath '\ori']);
    [~, maps.ori_ratio_polar] = GeneratePolarMap (maps.ori_ratio_params, maps.max_ori_ratio_change, 'colorshift', colorShift, 'outputpath', [outputpath '\ori']);
    [~, maps.ori_ratio_HLS_hc] = GenerateHLSMap(maps.ori_ratio_params, maps.max_ori_ratio_change, 'colorshift', colorShift, 'maxtuning', maps.tune_max, 'outputpath', [outputpath '\ori']);
    [~, maps.dir_ratio_HLS_hc] = GenerateHLSMap(maps.dir_ratio_params, maps.max_dir_ratio_change, 'colorshift', colorShift, 'maxtuning', maps.tune_max, 'outputpath', [outputpath '\dir']);
    imwrite(scaleImage(mean(data,3)),fullfile(outputpath,'avg.tif'), 'tiff');
else
    maps.dir_angle = GenerateAngleMap (maps.dir_dF_params.th, 'colorshift', colorShift);
    maps.ori_angle = GenerateAngleMap (maps.ori_dF_params.th, 'colorshift', colorShift);
    [~, maps.ori_ratio_polar] = GeneratePolarMap (maps.ori_ratio_params, maps.max_ori_ratio_change, 'colorshift', colorShift);
    [~, maps.ori_ratio_HLS_hc] = GenerateHLSMap(maps.ori_ratio_params, maps.max_ori_ratio_change, 'colorshift', colorShift, 'maxtuning', maps.tune_max);
    [~, maps.dir_ratio_HLS_hc] = GenerateHLSMap(maps.dir_ratio_params, maps.max_dir_ratio_change, 'colorshift', colorShift, 'maxtuning', maps.tune_max);
end



% dir_ratio_HLS = GenerateHLSMap(dir_ratio_params, max_dir_ratio_change, 'maxtuning', maps.tune_max);
% ori_dF_HLS = GenerateHLSMap (ori_dF_params, max_ori_dF, 'maxtuning', maps.tune_max, 'colorshift', colorshift, 'outputpath', [outputpath '\ori']);
% dir_dF_HLS = GenerateHLSMap (dir_dF_params, max_dir_dF, 'maxtuning', maps.tune_max, 'colorshift', colorshift, 'outputpath', [outputpath '\dir']);
% dir_dF_HLS_hc = GenerateHLSMap (dir_dF_params, max_dir_dF, 'maxtuning', maps.tune_max, 'HighContrast', 1);
% ori_ratio_HLS = GenerateHLSMap(ori_ratio_params, max_ori_ratio_change, 'maxtuning', maps.tune_max);


